{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "adjusted-humidity",
   "metadata": {},
   "source": [
    "## Validating Flash Calculations"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "manual-interpretation",
   "metadata": {},
   "source": [
    "Finding the solution to multiphase equilibrium calculations is challenging and the topic of continuing research.\n",
    "Many commercial packages offer users a great deal of confidence in their answers, but can they be trusted?\n",
    "Thermo can be used to validate the results from other software or identify defects in them."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "sitting-hygiene",
   "metadata": {},
   "source": [
    "The following example uses a natural gas mixture two pseudocomponents C7-C16 and C17+. The properties of pure components are taken from Thermo. To do a perfect comparison, the critical properties from other software packages should be substituted into Thermo.\n",
    "This is example S3 from Fonseca-Pérez (2021). The kijs are from Harding and Floudas (2000), and the original pseudocomponents are from Nagarajan, Cullick, and Griewank (1991).\n",
    "\n",
    "\n",
    "Fonseca-Pérez, R. M., A. Bonilla-Petriciolet, J. C. Tapia-Picazo, and J. E. Jaime-Leal. “A Reconsideration on the Resolution of Phase Stability Analysis Using Stochastic Global Optimization Methods: Proposal of a Reliable Set of Benchmark Problems.” Fluid Phase Equilibria 548 (November 15, 2021): 113180. https://doi.org/10.1016/j.fluid.2021.113180.\n",
    "\n",
    "Harding, S. T., and C. A. Floudas. “Phase Stability with Cubic Equations of State: Global Optimization Approach.” AIChE Journal 46, no. 7 (July 1, 2000): 1422–40. https://doi.org/10.1002/aic.690460715.\n",
    "\n",
    "Nagarajan, N. R., A. S. Cullick, and A. Griewank. “New Strategy for Phase Equilibrium and Critical Point Calculations by Thermodynamic Energy Analysis. Part I. Stability Analysis and Flash.” Fluid Phase Equilibria 62, no. 3 (January 1, 1991): 191–210. https://doi.org/10.1016/0378-3812(91)80010-S.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "french-mauritius",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "There are 2 phases present\n",
      "Mass densities of each liquid are 527.867861 and 527.867861 kg/m^3\n"
     ]
    }
   ],
   "source": [
    "from thermo import *\n",
    "from scipy.constants import atm\n",
    "pure_constants = ChemicalConstantsPackage.constants_from_IDs(\n",
    "    ['methane', 'ethane', 'propane', 'n-butane', 'n-pentane', 'n-hexane'])\n",
    "\n",
    "pseudos = ChemicalConstantsPackage(Tcs=[606.28,825.67], Pcs=[25.42*atm, 14.39*atm],\n",
    "                                   omegas=[0.4019, 0.7987], MWs=[140.0, 325.0])\n",
    "constants = pure_constants + pseudos\n",
    "\n",
    "properties = PropertyCorrelationsPackage(constants=constants)\n",
    "\n",
    "T = 353\n",
    "P = 38500e3\n",
    "zs = [0.7212, 0.09205, 0.04455, 0.03123, 0.01273, 0.01361, 0.07215, 0.01248]\n",
    "\n",
    "kijs = [[0.0, 0.002, 0.017, 0.015, 0.02, 0.039, 0.05, 0.09],\n",
    "[0.002, 0.0, 0.0, 0.025, 0.01, 0.056, 0.04, 0.055],\n",
    "[0.017, 0.0, 0.0, 0.0, 0.0, 0.0, 0.01, 0.01],\n",
    "[0.015, 0.025, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],\n",
    "[0.02, 0.01, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],\n",
    "[0.039, 0.056, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],\n",
    "[0.05, 0.04, 0.01, 0.0, 0.0, 0.0, 0.0, 0.0],\n",
    "[0.09, 0.055, 0.01, 0.0, 0.0, 0.0, 0.0, 0.0]]\n",
    "\n",
    "eos_kwargs = dict(Tcs=constants.Tcs, Pcs=constants.Pcs, omegas=constants.omegas, kijs=kijs)\n",
    "\n",
    "gas = CEOSGas(PRMIX, eos_kwargs, HeatCapacityGases=properties.HeatCapacityGases, T=T, P=P, zs=zs)\n",
    "liq = CEOSLiquid(PRMIX, eos_kwargs, HeatCapacityGases=properties.HeatCapacityGases, T=T, P=P, zs=zs)\n",
    "liq2 = CEOSLiquid(PRMIX, eos_kwargs, HeatCapacityGases=properties.HeatCapacityGases, T=T, P=P, zs=zs)\n",
    "phase_list = [gas, liq, liq]\n",
    "\n",
    "\n",
    "flashN = FlashVLN(constants, properties, liquids=[liq, liq2], gas=gas)\n",
    "# flashN.PT_SS_TOL = 1e-18\n",
    "res = flashN.flash(T=T, P=P, zs=zs)\n",
    "print('There are %s phases present' %(res.phase_count))\n",
    "print('Mass densities of each liquid are %f and %f kg/m^3' %(res.liquid0.rho_mass(), res.liquid0.rho_mass()))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "declared-appreciation",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The maximum relative difference in fugacity is 2.773985e-07.\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "max_fugacity_err = np.max(np.abs(1-np.array(res.liquid0.fugacities())/res.liquid1.fugacities()))\n",
    "print('The maximum relative difference in fugacity is %e.' %(max_fugacity_err,))"
   ]
  }
 ],
 "metadata": {
  "language_info": {
   "name": "python"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}